Setup

Loading and preparing data

### Import using haven
dat <-
  haven::read_sav(file.path(dataPath,
                            'database.orig.new labels.sav'));

### ... But convert to a data frame instead of a tibble, because,
### for example, length(unique(dat[, 'group'])) results in 1 instead
### of 2...
dat <- as.data.frame(dat);
dat <- haven::zap_label(dat);
dat <- haven::zap_labels(dat);

### Set group as factor
dat$group <-
  factor(dat$group,
         levels=1:2,
         labels=c("Alternative medicine",
                  "Medicine"));

dat$group_tri <-
  ifelse(dat$group=="Medicine",
         1,
         ifelse(dat$TCM6_simult < 3,
                2,
                ifelse(dat$TCM6_simult > 2,
                       3,
                       NA)));

dat$group_tri <-
  factor(dat$group_tri,
         levels=1:3,
         labels=c("Biomed",
                  "Alternative",
                  "Complementary"));
         
### Extract and store attitude variable names
attitudeVars <-
  grep('ATT1_',
       names(dat),
       value=TRUE);

attitudeVars_2 <-
  grep('ATT2_',
       names(dat),
       value=TRUE);

Construct tree of attitude items

constructTreeYAML <-
  yum::load_yaml_fragments(here::here("methods-construct-tree",
                                      "cam-biomed-attitude-tree-1.dct"));

constructTree <-
  yum::build_tree(constructTreeYAML);

constructTree$Do(function(node) {
  nameFromDataset <-
    grep(node$name,
         names(dat),
         value=TRUE);
  print(nameFromDataset);
  if (length(nameFromDataset) > 0) {
    node$label <- nameFromDataset;
  }},
  filterFun = data.tree::isLeaf);

## [1] "ATT1_4_BodyMirrorToSoul"
## [1] "ATT1_31_IllnessSymbology"
## [1] "ATT1_21_IllnessBecauseOfEmotions"
## [1] "ATT1_33_SymptomsWillDisappear"
## [1] "ATT1_9_MindHasStrongEffect"
## [1] "ATT1_11_BodyRemembers"
## [1] "ATT1_36_UnprocessedTrauma"
## [1] "ATT1_28_ConfrontingEmotionalProb"
## [1] "ATT1_26_HealingResultOfEmoDevelopment"
## [1] "ATT1_18_TreatsOnlySymptoms"
## [1] "ATT1_17_SickByChance"
## [1] "ATT1_20_HealingIsLuck"
## [1] "ATT1_10_AttractPplAndEvents"
## [1] "ATT1_16_EverythingConnected"
## [1] "ATT1_23_NothingByChance"
## [1] "ATT1_2_MustSuffer"
## [1] "ATT1_38_IllnessTeachesUs"
## [1] "ATT1_34_EnergyInEastern"
## [1] "ATT1_27_EnergeticSystemInBody"
## [1] "ATT1_14_Reincarnation"
## [1] "ATT1_6_IllnessIsImbalance"
## [1] "ATT1_40_StrongerComplaintsMeanHealing"
## [1] "ATT1_29_RadiationTherapyHarmful"
## [1] "ATT1_25_ChemotherapyHarmful"
## [1] "ATT1_39_OnlyNatural"
## [1] "ATT1_22_AvoidPharma"
## [1] "ATT1_37_NoBiopsy"
## [1] "ATT1_30_MandatoryVaccines"
## [1] "ATT1_15_TrustAncientRemedies"
## [1] "ATT1_3_TrustInTradRemedy"
## [1] "ATT1_13_ElectronicRadiation"
## [1] "ATT1_7_HealthyDiet"
## [1] "ATT1_1_ExerciseAndDiet"
## [1] "ATT1_12_WeakImmuneSystem"
## [1] "ATT1_24_IllnessBecauseOfGenes"
## [1] "ATT1_32_TrustWesternDocs"
## [1] "ATT1_35_SeriousSymptom"
## [1] "ATT1_5_DependsOnEnvironment"
## [1] "ATT1_19_DoctorMustHealMe"
## [1] "ATT1_8_NeedTestResult"

### Set labels as names
constructTree$Do(function(node) node$name <-
                   node$label);

### Convert to DiagrammeR graph
constructGraph <-
  data.tree::ToDiagrammeRGraph(constructTree);

### Show graph
DiagrammeR::render_graph(constructGraph);

### Export graph
DiagrammeR::export_graph(constructGraph,
                         file_name = here::here("methods-construct-tree",
                                                "cam-biomed-attitude-tree-1.png"));

### Also plot as dendrogram (method not exported by this version of data.tree, oddly)
constructDendro <-
  data.tree:::as.dendrogram.Node(constructTree);

### Get labels in same order
constructTreeLabels <-
  unlist(constructTree$Get('label', filterFun=data.tree::isLeaf));

### For future reference: check
### http://www.sthda.com/english/wiki/beautiful-dendrogram-visualizations-in-r-5-must-known-methods-unsupervised-machine-learning#ggdendro-package-ggplot2-and-dendrogram

ggConstructDendro1 <-
  ggdendro::ggdendrogram(constructDendro,
                         rotate=TRUE,
                         theme_dendro = TRUE) +
  ggplot2::scale_x_continuous(position="top",
                              breaks=seq_along(constructTreeLabels),
                              labels=constructTreeLabels) +
  ggplot2::scale_y_reverse();

## Scale for 'x' is already present. Adding another scale for 'x', which
## will replace the existing scale.

## Scale for 'y' is already present. Adding another scale for 'y', which
## will replace the existing scale.

print(ggConstructDendro1);

ggsave(filename=here::here("methods-construct-tree",
                           "cam-biomed-dendrogram-1.png"),
       plot=ggConstructDendro1,
       width=12,
       height=19,
       units='cm');

ggConstructDendro2 <-
  constructDendro %>%
    dendextend::set("branches_k_color",
                    value = viridis::viridis(4),
                    k = 4) %>%
    dendextend::as.ggdend() %>%
    ggplot2::ggplot(horiz=TRUE);

print(ggConstructDendro2);

## Warning: Removed 69 rows containing missing values (geom_point).

ggsave(filename=here::here("methods-construct-tree",
                           "cam-biomed-dendrogram-2.png"),
       plot=ggConstructDendro2,
       width=40,
       height=30,
       units='cm');

## Warning: Removed 69 rows containing missing values (geom_point).

Descriptives

ufs::cat0("\n\n### Group\n\n");

pander(userfriendlyscience::freq(dat$group));

ufs::cat0("\n\n### Missing values\n\n");

apply(is.na(dat[, attitudeVars]), 2, sum);

           ATT1_1_ExerciseAndDiet 
                                0 
                ATT1_2_MustSuffer 
                                0 
         ATT1_3_TrustInTradRemedy 
                                0 
          ATT1_4_BodyMirrorToSoul 
                                0 
      ATT1_5_DependsOnEnvironment 
                                0 
        ATT1_6_IllnessIsImbalance 
                                0 
               ATT1_7_HealthyDiet 
                                0 
            ATT1_8_NeedTestResult 
                                0 
       ATT1_9_MindHasStrongEffect 
                                0 
      ATT1_10_AttractPplAndEvents 
                                0 
            ATT1_11_BodyRemembers 
                                0 
         ATT1_12_WeakImmuneSystem 
                                0 
      ATT1_13_ElectronicRadiation 
                                0 
            ATT1_14_Reincarnation 
                                0 
     ATT1_15_TrustAncientRemedies 
                                0 
      ATT1_16_EverythingConnected 
                                0 
             ATT1_17_SickByChance 
                                0 
       ATT1_18_TreatsOnlySymptoms 
                                0 
         ATT1_19_DoctorMustHealMe 
                                0 
            ATT1_20_HealingIsLuck 
                                0 
 ATT1_21_IllnessBecauseOfEmotions 
                                0 
              ATT1_22_AvoidPharma 
                                0 
          ATT1_23_NothingByChance 
                                0 
    ATT1_24_IllnessBecauseOfGenes 
                                0 
      ATT1_25_ChemotherapyHarmful 
                                0 

ufs::cat0("\n\n### Attitude\n\n");

ufs::meansComparisonDiamondPlot(dat,
                                rev(attitudeVars),
                                compareBy = 'group_tri',
                                comparisonColors = viridis::viridis(3,
                                                                    end=.7),
                                dataAlpha=.25);

Network analyses

For the network analysis section we will first estimate a full network where a choice of medical practice is regarded as a system component. This will allow us to explore how attitudes relate to the behavior of interest and the relative importance of each determinant. Afterwards, we will create sub samples based on a variable that represents groups of people that prefere one or other medical practice. This way we can investigate structural features of networks (topologies) that are peculiar to these groups.

# here I subset data that is going to be used in the network analysis with the behavior as a component.

subset <- dat[, c(attitudeVars, 
                 "group_tri")];

# I eliminated other variables of interest (sex, age etc.) because the N in subgroups are very small and the number of "predictor" variables would exceed the number of datapoints.  


### Here I create subsets based on the grouping variable.
for(i in levels(subset$group_tri)){
  assign(paste("subset",
               i,
               sep = "_"),
         subset(subset,
                subset$group_tri == i))    
}

subset$group_tri <-
  as.numeric(subset$group_tri);

subset <-
  subset %>% tidyr::drop_na();

Step 1: Estimating the netork with full data.

network <-
  estimateNetwork(subset,
                  default = "EBICglasso");

## Estimating Network. Using package::function:
##   - qgraph::EBICglasso for EBIC model selection
##     - using glasso::glasso
##   - qgraph::cor_auto for correlation computation
##     - using lavaan::lavCor

## Variables detected as ordinal: ATT1_1_ExerciseAndDiet; ATT1_2_MustSuffer; ATT1_3_TrustInTradRemedy; ATT1_4_BodyMirrorToSoul; ATT1_5_DependsOnEnvironment; ATT1_6_IllnessIsImbalance; ATT1_7_HealthyDiet; ATT1_8_NeedTestResult; ATT1_9_MindHasStrongEffect; ATT1_10_AttractPplAndEvents; ATT1_11_BodyRemembers; ATT1_12_WeakImmuneSystem; ATT1_13_ElectronicRadiation; ATT1_14_Reincarnation; ATT1_15_TrustAncientRemedies; ATT1_16_EverythingConnected; ATT1_17_SickByChance; ATT1_18_TreatsOnlySymptoms; ATT1_19_DoctorMustHealMe; ATT1_20_HealingIsLuck; ATT1_21_IllnessBecauseOfEmotions; ATT1_22_AvoidPharma; ATT1_23_NothingByChance; ATT1_24_IllnessBecauseOfGenes; ATT1_25_ChemotherapyHarmful; ATT1_26_HealingResultOfEmoDevelopment; ATT1_27_EnergeticSystemInBody; ATT1_28_ConfrontingEmotionalProb; ATT1_29_RadiationTherapyHarmful; ATT1_30_MandatoryVaccines; ATT1_31_IllnessSymbology; ATT1_32_TrustWesternDocs; ATT1_33_SymptomsWillDisappear; ATT1_34_EnergyInEastern; ATT1_35_SeriousSymptom; ATT1_36_UnprocessedTrauma; ATT1_37_NoBiopsy; ATT1_38_IllnessTeachesUs; ATT1_39_OnlyNatural; ATT1_40_StrongerComplaintsMeanHealing; group_tri

plot(network,
     layout = 'spring',
     labels = colnames(network),
     title = c('Figure 1: A Network with the behavioral variable'));

### Also store to disk - Sam, you can copy this for any other plots of course
pdf(here::here('results-intermediate',
               'network-with-grouping-variable.pdf'));
plot(network,
     layout = 'spring',
     labels = colnames(network),
     title = c('Figure 1: A Network with the behavioral variable'));
dev.off();

png 2

# Calculating centrality measures.
pdf(here::here('results-intermediate',
               'centrality_plot.pdf'));
centralityPlot(network,
               include = "all");

## Note: z-scores are shown on x-axis rather than raw centrality indices.

dev.off();

png 2

# # Checking the stability of the centrality measures
# central_stability <- 
#   bootnet(network,
#           nCores = 20,
#           nBoots = 1000,
#           type = 'case');
# 
# pdf(here::here('results-intermediate',
#                'centrality_stability_plot.pdf'));
# plot(central_stability)
# dev.off();
# 
# # Checking the stability/reliability of the edge weights
# edgewgt <-
#   bootnet(network,
#           nCores = 20,
#           nBoots = 2500);
# 
# plot(edgewgt,
#      labels = FALSE,
#      order = 'sample');
# 
# pdf(here::here('results-intermediate',
#                'edge_weights.pdf'));
# plot(edgewgt,
#      labels = FALSE,
#      order = 'sample');
# dev.off();

### First we need to take the absoulte values of the adjacency matrix and devide 1 by the matrix. Then recreate a network object for further analysis.

absolute_adj <-
  abs(network$graph);

for (i in which(absolute_adj > 0)) {
  absolute_adj[i] = 1/absolute_adj[i] 
}

graph_full <- 
  graph.adjacency(absolute_adj,
                  mode = 'undirected',
                  weighted = TRUE);

# Calculate shortest path to the outcome variable and then delete the last row of the dataframe that includes the outcome variable (shortest path to itself = 0) 

dijkstra_fullnetwork <-
  igraph::distances(graph_full,
                    v = V(graph_full),
                    to = 41,
                    algorithm = "dijkstra");

dijkstra_ful <-
  subset(dijkstra_fullnetwork,
         dijkstra_fullnetwork[,] == min(dijkstra_fullnetwork[-c(41),]));

fg <- 
  fastgreedy.community(graph_full, weights = E(graph_full)$weight)

fg$names <- 
  strtrim(fg$names, 7)

V(graph_full)$name <-
  strtrim(fg$names, 7)

length(fg)

sizes(fg)

set.seed(100)

par(mfrow=c(1,2))

plot(fg, graph_full, 
     vertex.label.cex=c(0.5,0.5,0.5),
     vertex.label.font=c(2))

dendPlot(fg, mode = 'phylo')

Step 2: Group Comparisons.

alternative <-
  tidyr::drop_na(subset_Alternative[,-c(41)])
biomedical <-
  tidyr::drop_na(subset_Biomed[,-c(41)])

network_fgl <-
  EstimateGroupNetwork(list('Alternative' = alternative,
                            "Biomed" = biomedical)) # We get an empty network.

absolute_adj_hc <-
  abs(network$graph);


graph_full_hc <- 
  graph.adjacency(absolute_adj[-c(41),-c(41)],
                  mode = 'undirected',
                  weighted = TRUE)


fg_item <- 
  fastgreedy.community(graph_full_hc, weights = E(graph_full_hc)$weight)

fg_item$names <- 
  strtrim(fg_item$names, 7)

V(graph_full_hc)$name <-
  strtrim(fg_item$names, 7)

length(fg_item)

sizes(fg_item)

set.seed(101)

par(mfrow=c(1,2))

plot(fg_item, graph_full_hc, 
     vertex.label.cex=c(0.5,0.5,0.5),
     vertex.label.font=c(2))

dendPlot(fg_item, mode = 'phylo')